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ABSTRACT 

We study how the magnetorotational instabihty (MRI) in protoplanetary disks is affected 
by the electric discharge caused by the electric field in the resistive MHD. We have performed 
three-dimensional shearing box simulations with various values of plasma beta and electrical 
breakdown models. We find the self-sustainment of the MRI in spite of the high resistivity. 
The instability gives rise to the large electric field that causes the electrical breakdown, and the 
breakdown maintains the high ionization degree required for the instability. The condition 
for this self-sustained MRI is set by the balance between the energy supply from the shearing 
motion and the energy consumed by the Ohmic dissipation. We apply the condition to various 
disk models and study where the active, self-sustained, and dead zones of MRI are. In the fiducial 
minimum-mass solar nebula (MMSN) model, the newly-found sustained zone occupies only the 
limited volume of the disk. In the late-phase gas-depleted disk models, however, the sustained 
zone occupies larger volume of the disk. 



Subject headings: Dust - 
— MHD — instabilities 

Introduction 



Protoplanetary disks are the sites of planet for- 
mation. The disk turbulence greatly affects the 
mutual sticking of the planetcsimals, their settle- 
ment to the disk midplanc. The turbulence is the 
source of the angular momentum transfer in the 
disk that causes gas accretion and migration of 
the planetesimal onto the central star. Thus un- 
derstanding the evolution of the turbulence within 
protoplanetary disks is an essential step both in 
the studies of the disk evolution and the planet 
formation. 



planets and satellites:formation — planetary systems: protoplanetary disks 



The magnetorotational instability (MRI) is 
considered to be the major source of turbulence 
in many types of accretion disks inc l uding pro- 
toplanetary disks (jBalbus fc HawlevI (|l998l ) and 
references therein). One of the distinct proper- 
tics of the protoplanetary disks compared to other 
accretion disks is that the major parts of the pro- 
toplanetary disks are only weakly ioniz ed, and the 



magne t ic diffusivity a f fects the MRI (jSano et al 



(| 19981) : [Fleming et al.l (|2000[) .) The low ionization 
degree is due to their low temperature and high 
number density of the dust component. 
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In protoplanetary disks, MRI and the dust com- 
ponents affect each other. The turbulence is one 
of the source for the relative veloci t ies of the col- 
liding dust (jOrmel fc Cuzzil 120071 : iBrauer et al. 
20081 ). and co ntributes both on dust growth 



and disruption dBlum fc Wurmll200l lWada et al 



20081 : iGiittler et al.l 120101: IWettlaufed l2010[ ). On 
the other hand, dust particles in protoplane- 
tary disks are the major sites of charged particle 
recombination, and thereby influences the ion- 
ization degree of the disk. dSodha et al. 2009t 
Umebavashi fc NakanolboOOl : ICrach et al.ll2010[ ) 

The dead zone can occupy a large volume of 
a protoplanetary disk, especially in the presence 
of abundant small dust grains (iGammid Il996t 
Sano et al.|[2000t lllgner fc ]N"elsonll2006l ). However, 
various electric discharge mechani sms in proto- 



plane t ary disks have be e n pro p osed (IHoranyi et al. 
199,4 iDesch fc Cu^ I2OOOI: iMuranushil I2OIOI ) 



which may provide higher ionization degree com- 
pared to the values predicted by the dust-absorption 
equilibrium models, resulting in the increased MRI 
activity in the disk. They consider the electron 
avalanche process, an exponential growth in the 
number of conducting electrons that takes place 
when the kinetic energy of the electrons exceeds 
the ionization energy of a neutral gas molecule. 
The result is electrostatic breakdown, the low- 
ering of the resistivity of the fluid and electric 
discharge, increase of the electric current through 
the fluid. 

Moreover, a model is proposed where the MRI 
itself provides sufficient ionization (jlnutsuka fc Sano 
2OO5I hereafter IS05). IS05 have shown that that 
the electric field typically generated by the proto- 
planetary disk turbulence is strong enough to drive 
the electrons away from the thermal Maxwell- 
Boltzmann distribution. Those energetic electrons 
contained in electric current cause the electric dis- 
charge and maintains the ionization degree high 
enough for the MRI to survive. IS05 have also 
shown that the energy supply from the shearing 
motion is about 30 times larger than the energy 
required to maintain the sufficient number of elec- 
trons in the presence of standard dust grains. 

However, IS05 have studied only one-zone mod- 
els, and with only one set of parameters typical to 
1 AU of the disk. In this work, we extend the 
model IS05 to a local, 3D simulations of proto- 
planetary disks and study the interaction of the 



MRI with the discharge ionization. We also ap- 
ply the model to global models of protoplanetary 
disks and study where and when in the disk the 
self-sustainment of MRI takes place. 

This paper is organized as follows, in §2, we 
perform the numerical simulations of the MRI 
in unstratified, th ree-dimensiona l shea ring-boxes 



along the lines of Hawlev et al. ( 19951 hereafter 



HGB95) with the nonlinear Ohmic diffusivity 
added. In §3, We analyze the activity of the MRI 
in the protoplane tary disk, using the method of 
Sano et al.l (|200d hereafter SMUNOO). To calcu- 
late the ionization de gree in the disk we use the 
method proposed by lOkuzumil (|2009l hereafter 
009). §4 is devoted to conclusions and discus- 
sions. Table [T] lists the symbols frequently used in 
this paper. 

2. Simulations of The MRI with nonlinear 
Ohm's law 

2.1. Numerical Setup 

There are three diffusion terms in MHD; they 
are Ohmic diffusion. Hall diffusion and ambipolar 
diffusion. In protoplanetary disks any one of the 
three modes can be the domin ant mode depe nding 
on dust and gas density (e.g. IWardld 2007 ). and 
interactio n between the different m odes may alter 
the MRI (jWardle fc Salmeronl20"li ) . In this paper 
we only focus on the Ohmic diffusion because it is 
the most studied one in the context of the MRI. 
We leave the treatment of other diffusion modes 
for future studies. 

The electric discharge taken into account, we 
construct a simple model of the discharge as fol- 
lows, in terms of an appropriate rjo and Jcrit^ 



E' = 

v{J) = 



47r 



J. 



crit 



J 



-'70 



if 
if 



J < Jr. 



J> J, 



crit ■ 



(1) 



(2) 



This nonlinear diffusivity model states that the 
electric field on the fluid co-moving frame never 
exceeds a critical value, £"crit = 47rc~^r7o Jd-it, 
thus the magnetic diffusivity rj varies depending 
on electric current J, and Ohm's law become non- 
linear. 

Note that the smallest space scale dealt in this 
paper is of order of 10~^ AU. The actual scale 
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Symbol 


Value (Dimension) 


Definition 


Location 


P 


(g cm ^) 




Gas density 


m 


V 


(cm s^-'-) 




Gas velocity 


m 


B 


(gi/2 cm~^/^ s" 


') 


Magnetic field 


m 


P 






Pressure with isothermal equation of state 


m 


E 


(gl/2 


') 


Electric field in the lab frame 


m 


E' 


(gl/2 


') 


Electric field in the comoving frame 


(??) 


J 


/ 1 /2 — ^/2 - 

[g ' cm ■^'^ s 


1 \ 

) 


Electric current 


COD 


Tin 

'/O 


(cm2 s-i) 








J crit 


(gi/2 cm~^/^ s~ 




— critical current for extended Ohm's law 


m 




(g cm~-^s^^) 




Initial pressure 




/3 






Plasma beta 










Initial, vertical net magnetic field 


(USD 


H 






Disk scale-height 


m 








The nondimensionalization unit of magnetic field 


m 








The nondimensionalization unit of current 


m 


Rm 






Magnetic Reynolds number 


mi 


h 


1.0 




Surface Density Multiplier 




q 


3/2 




Power Law Index of the Surface Density 




ad 


0.1 /i m 




Radius of Solid Dust Particle 


(gSl) 


fd 


0.01 




Dust to Gas Ratio 





Table 1: The list of symbols used in this paper. 



of the discharge structures can be much smaller 
than this. The estimate for the lower limit of the 
size of such structures is their thickness, which is 
of the order cif 500 times electron mean free path 
(|PilipD et al.lll992l ). 

However, we can derive the macroscopic dis- 
charge model from the microscopic discharge re- 
lation \E'\ < E'crit that holds everywhere in the 
plasma. Since e.g. the x-component of the dis- 
cretized electric field {E'^) is obtained from the 
line integral of the real field over the discretiza- 
tion length AV, 



\{E'.)\ 



J E'.dV 



< 



< 



AV 

I\E'.\dV 

AV 
J E'„i,dV 

AV 

E crit- 



(3) 



Therefore, we can use Eqs. ([T]), ([2]) as a "coarse- 
grained model" where we can interpret E', J as 
spatial averages. If electrical breakdown occurs in 



a scale smaller than grid size, the spatially aver- 
aged electric field is smaller than £"crit- Thus, in 
general, the electrical breakdown may occur even 
in the region where (£"crit) is smaller than £"crit- 
Therefore, if we adopt equations ([T|) and ((2)), we 
may underestimate, but not overestimate, the oc- 
currence of electric discharges. 

We adopted a local, Kepler-rotation shearing 
box that has radial (x), azimuthal (y) and vertical 
(z) axes, and solved the following resistive magne- 
tohydrodynamic (MHD) equations numerically: 



|+V.(pv) = 0, 
dv 



dt 



V • Vv 



1 



-V P 



5B 

'dt ^ ' 
with isothermal EOS 

P = 



P \ Stt, 
cV X E, 



cIp,. 



47rp 



(4) 

(B • VB) 

(5) 
(6) 



(7) 
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Had we used the adiabatic EOS, the internal en- 
ergy would have kept growing as the linear func- 
tion of time (HGB95). Therefore, we use the 
isothermal EOS to approximate the steady state 
attained by the cooling processes present in the 
protoplanetary disks. 

The nonlinear Ohm's law reads: 

E = -ivxB + ^?/(J)J, (8) 
c & 

(9) 
(10) 



J = —V X B, 

47r 



Here, we have studied three different models 
for nonlinear diffusivity r/(J). In addition to our 
fiducial model(f id), Eqs. ^ ([2]), we have studied 
the following two models: 



(p2): ,7 (J) 



(p4): ,7 (J) = 1 



J crit 
J crit 



,(11) 



-4\ 4 



(12) 



The electric field as functions of current density in 
these three models are shown in Figure [T] 

Following HGB95, we set up our numerical ini- 
tial conditions as follows. We use the disk scale- 
height H as the unit length. The box size is 
{Lx, Ly, Lz) = (l,27r, 1). First, we set the aver- 
age values po = 1 and Pq = 10^^ to every mesh, 
and let fluid velocity to be at rest in shearing 
box frame; {vx,Vy,Vz = 0, — (3/2)rix, 0). Here, 
Cs = 10"^ and also Vt = lO^^. 

Next, we introduce random perturbations in 
density, pressure, and velocity. The density and 
pressure perturbations are in proportion so that 
the isothermal EOS is met, and the amplitude 
is 5p/po = SP/Po = 2.5 X 10-2. We perturb 
the velocity component-wise, with the amplitude 
Svi = 5 X 10~^Cs for each. 

We use the following units of magnetic field and 
electric current: 



and the scale height 
H 



rR 

AttH ' 



n' 



(13) 
(14) 

(15) 



We set uniform magnetic field in the z-direction, 
and express the initial field strength by the plasma 
beta. 



/3 = SeqpVS.o' - 87rPo/S,o'. 



(16) 



The plasma beta satisfy the following relation 
between the sound speed Cs and the Alfven veloc- 
ity along the magnetic field • 



/3 



2c2 



(17) 



We define the magnetic Reynolds number as: 

Rm = v^^y,^o^, (18) 



using the Alfven velocity v^^ = S^o/V^Trp set 
by the initial vertical magnetic field. This is in 
accordance with SMUNOO and IS05 while some 
literature a dopts different defiii ition (e.g. Rm = 
Cs^/rjofl in Fleming et al. 2000[ ) . 

We have used Athena ( Gardiner fc Ston3 2005 . 



2008[ ). an open-source MHD code for our simula- 
tions. 

2.2. Simulations Procedure 

We vary the initial magnetic field strength and 
the diffusivity models, and we classify each set 
of parameters as either O^-ctive zone, xdead 
zone, or Asustained zone. The experiment 
method and the definition of the three classes are 
given in this section. 

The parameters we have investigated are the 
initial vertical field strength (represented by 
plasma /?), the linear diffusivity (represented by 
magnetic Reynolds number Rm), and the criti- 
cal current Jdit- The range of the survey was 
400 < ^ < 25600, 0.002 < i?M < 2 and 
0.01 < Jcrit/>/cqp < 100. In addition, the lim- 
iting cases of Rm = oo and Jcnt/ Jcqp = oo, that 
respectively correspond to ideal MHD models and 
linear Ohm's law models, are studied for compar- 
ison with the literature. 

For each value of (3, we prepared the initial con- 
dition as described in section 12.11 and continued 
the simulation for 10 orbits (t ~ 207r/r2), at first 
with magnetic diffusivity turned off {ri(J) = 0). 
While running the simulations, we created restart 
data for every periodic points {t = inn /Q. where n 



4 



is an integer). The condition allowed the MRI to 
grow and saturate in about 5 orbits {t = IOtt/JI). 

Then, for each pair of (i?M, Jait/ Jcqp), we 
turned on the difFusivity and re-started the sim- 
ulation either from the initial laminar flow (t ~ 
0) or the saturated MRI states at 8, 9 and 10 
orbit {t = 167r/f}, 187r/17,207r/17). We numeri- 
cally evolved them until they reach 20 orbit (t = 
407r/r2). The reason why we have adopted three 
different MRI saturated initial conditions (t = 
167r/ri, 187r/r2, 207r/f2) for each set of the param- 
eters is that a 'turbulent initial condition' is not 
unique; therefore we need to test if our results de- 
pend on the choice of the initial condition or not. 

During each simulation run, we recorded the 
space averages of physical quantities as the func- 
tions of time, such as magnetic energy density 
B^, the Reynolds and Maxwell stress pv^-Svy, 
—BxBy/4:TT, and the squared current J^. After 
the simulations we studied the time average of the 
quantities. For a physical quantity A, we denote 
its space and time average by (A) and A, respec- 
tively. Their definitions are as follows: 



(A) 
A 



J dx J dy J dzA 
Jdxjdyjdz 
JdtA 
-Jdf- 



(19) 
(20) 



The space average is taken for the entire compu- 
tational domain (0 < x < L^., < y < Ly, < 
z < Lz): numerical resolution is {N^, Ny, Nz) = 
(64, 128, 64), and the time average is taken for the 
last five orbit (SOtt < t < 407r) unless otherwise 
mentioned. The statistics for the important sets 
of parameters are presented at the end of the pa- 
per. 

Using the average values, we classify each set 
of the parameter (/?, Rm, ^crit/^eqp) as follows 
(c.f. Table [2|). First, a parameter is in Q)a.ctive 
zone if the MRI is observed both in the simula- 
tion started from laminar flow as well as in all of 



zone 


from laminar? 


from ideal MRI? 


Qactive 


unstable 


unstable 


xdcad 


stable 


stable 


A sustained 


stable 


unstable 



Table 2: The zone names and the meaning of the 
symbols in Figures |3] and 



the three simulations started from the saturated 
MRI states. Second, a parameter is in xdead 
zone if the instability is not observed neither in 
the simulation started from laminar flow as well 
as in any of the three simulations started from the 
saturated MRI states. Finally, a parameter is in 
Asustained zone, if the MRI is observed at in all 
the three simulations started from the saturated 
MRI states, but not at the end of the simulation 
started from laminar flow. 

2.3. The Result of Shearing-Box Simula- 
tions 

The typical behavior of the current for the ac- 
tive zone, dead zone and sustained zone are in Fig- 
ure [H From the simulations we have observed that 
the three classes (active, dead and sustained) are 
exhaustive: that the runs started from 8,9 and 10 
orbit always agree in terms of the classification; 
and that if the MRI dies when starting from sat- 
urated initial condition, it also does not activate 
starting from laminar initial condition. 

To classify the active, dead and sustained zone, 
we need to assess the magnctorotational instabil- 
ity of the system, so we introduce the following 
criteria for quantitative assessment. We say that 
the system is magnetorotationally unstable if the 
averaged current (J^)^/^ is greater than O.lJoqp, 
and stable if otherwise. Here, the time average 
is taken for the last five orbit (SOtt < t < AOtt). 
We have learned from the simulations that the 
quantity (J^)^/^ is a good indicator for the sta- 
bility, since it either fluctuates around mean value 
(j2)i/2 ^ loj^qp (unstable) or go under (j2)i/2 < 
0.1 Jeqp almost monotonically with little vibration 
(stable), and there is no ambiguity between the 
two; c.f. Figure [2] 

However, as we make Rm smaller, the diffu- 
sion timcscale becomes shorter, and the wall clock 
time for the simulations until t = 407r becomes 
impractically larger. Therefore, for the parameter 
range Rm < 0.01, we terminate the simulations at 
1.5 X 10^ cycles and determine the class by extrap- 
olations. 

Figures [3] and [5] show the result of our param- 
eter survey. We found that sustained zone exist 
— the MRI does exhibit hysteresis behavior for a 
certain set of parameters. 

To study possible influences of the numeri- 
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Fig. 2. — Time evolution of the averaged current density for three magnetic diffusivity models. Models are 
(a) /3 = 400, Rm = 0.6, Jcrit/Jcqp = 1 (b) /3 = 400, Rm = 0.2, Jcrit/Jcqp = 1 (c) /3 = 400, Rm = 0.2, 
Jcnt/Jcqp = 10 . The graphs show typical current behavior for (a) Q active zone, (b) A sustained zone, or 
(c) X dead zone. 
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Fig. 3. — Distribution of the sets of parameters (/3, Rm, Jmt/ Jccip) in our simulations. We classify each 
of them as either Q active zone, A sustained zone, or x dead zone. This page includes the data for 
400 < /3 < 3200. The parameters classified by extrapolations are marked by light gray symbols. 
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Fig. 4.— Continued from Figure El the data for 6400 < /3 < 25600. 
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Fig. 5. — Dependence of the MRI behavior on the rcsohition. The development of electric current over time 
in the simulations using the different numerical resolution Nx ~ 32, 48, 64(fiducial) and 80 while keeping the 
aspect ratio : Ny : iV^ = 1 : 2 : 1, for three different set of parameters (a) (3 = 400, Rm = 0.6, 
Jcrit/Jeqp = 1 , (b) /3 = 400, Rm = 0.2, Jcrit/Jcqp = 1 and (c) /3 = 400, Rm = 0.2, Jcrit/Jcqp = 10 . They 
are typical parameter sets for (a) active zone, (b) A sustained zone, and (c) x dead zone, respectively. 
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Fig. 6. — Dependence of the MRI behavior on resistivity models. The development of electric current over 
time in the simulations using different resistivity models f id( Eqs. ([T]) ([2])), p2 (Eq. [TT|) and p4 (Eq. [T2|l . 
for three sets of parameters (a) /3 = 400, Rm = 0.6, JcHt/Jcqp = 1 j (b) /3 = 400, Rm = 0.2, Jcnt/Je^p = 1 
and (c) /3 = 400, Rm = 0.2, Jciit/Joqp = 10 ; which are typical parameter sets for (a) Q active zone, (b) A 
sustained zone, and (c) x dead zone, respectively. 
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cal resolutions, we have performed simulations 
using the different numerical resolution Nx = 
32, 48, 64(fLducial) and 80 while keeping the as- 
pect ratio Nx : Ny : TV^ = 1 : 2 : 1, for three 
different set of parameters {P, Rm , Jcrit/ Jcqp) = 
(400,0.6,1), (400,0.2,1) and (400,0.2,1). Figure [5] 
summarizes the convergence tests. Within these 
sets of parameters, we observe that our classifica- 
tion does not depend on the numerical resolution. 

We have also studied if our result depends on 
the model of nonlinear Ohm's law. In addition to 
our fiducial(f id) model, Eq. ([2]), we have studied 
two smoothly transitting models, Eqs. (fTTj) and 
([T2|) . To distinguish the three models see Figure 



tion (EU): 



m 

Figure [6] summarizes the time evolution of the 
current density in the simulations with the three 
typical sets of parameters {P,Rm, J mt/ J cqp) = 
(400,0.6,1), (400,0.2,1) and (400,0.2,1). Within 
the regime we have tested, the hysteresis behav- 
ior does not depend on the detail of the non-linear 
resistivity models. 

From Figures [3] and [4] we can see the following 
condition for the Asustained zone: 



J. 



crit 



1 



hb, 



(21) 



where /whb is a proportionality constant that sat- 
isfies /whb 5- 15 for 400 < /3 < 1600, and 
/whb — 15 — 50 for /3 > 3200. Hereafter, we inter- 
pret this in terms of the work-heat balance in the 
resistive MRI. 

We also remark that within parameter re- 
gions that is in active zone, Rm < 1, and with 
larger Jcrit, we observed the large- amplitude time- 
variability of physical quantities such as magnetic 
fields due to repeated growth and reconncction 
of the cha nnel solutions. This phenomenon is re- 
ported by [Fleming et al.l (|2000l) . 



2.4. 



Interpretation of The Simulation Re- 
sults 



In this section, we show that Equation ([?!]) 
can be understood as a condition of the balance 
between the magnetic energy dissipated by Joule 
heating per unit volume (Wj) and the work done 
by shearing motion per unit volume (Wsh)- 

Let us define Qwhb as the left hand side of Equa- 



Q 



whb 



•/crit 1 
Jcnn Rm 



(22) 



The condition for self-sustained MRI is Qwhb < 
/whb; which is to be explained in this section. 



First, substitute Rm = ^^Az^/'yo^; 



i whb 



•/crit ^70^^ 
7 V 2 



(23) 



Next, when the MRI is active, (J^)^/^ ~ 
./sat'/cqp where /sat is of the order of 10. Since 
this average current strength lies in super-critical 
regime of Ohm's law (J > Jcrit): the electric field 
is E' crit — 47rc~^77o'^rit as modeled in Equation 
dl]) and dH). Therefore, Joule heating per unit 
volume is estimated as 



— E clit ■ fsaxJcqp 



47r/s, 



"Ho J crit J 



(24) 



To explain the work-heat balance qualitatively, 
this estimate needs correction due to the high 
space variability of current field under the dis- 
charge conditions. We introduce /fin, the filling 
factor, the ratio of the volume that contributes to 
the Joule heating to the total volume. Formally, 
/fill is defined as the ratio between volume aver- 
ages of the actual Joule heat generated and the 
Joule heat estimated by this method: 



/fill = 



{E'J) 



(j2)l/2 



(/(■/) J) 
(j2)l/2 ■ 



where 



nj) 



if 
if 



J > t/crit 
J < Jcrit 



(25) 



(26) 



Using this /hn. Equation is rewritten as: 

Wj = 47r/fiii/satC"^?7oJcritJcqp- (27) 

Substituting 770 in with ^7} gives 



ij whb 



w.i 



47r/fill/satJoqp^WAz^ ' 

47r /fill /sat >/cqp ^ ^AZ ^ Q whb 



(28) 
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Substituting v^^ with Equation (|17p and then 



Cs with ([T5|) . Jcqp with Equation (fT4|) and then 
-Beqp with Equation (fT3| gives 



whb 



On the other hand, 



3 



(29) 



(30) 



where a is lShakura fc SunvaevI (|1973[ ) 's a param- 
eter. Substituting Equation (|30|) into Equation 
(l29l). one obtains 



8 /fill /s 



sat ^; whb 



3a/3 



Vt^sh. 



(31) 



For the MRI to sustain itself by the discharge pro- 
cess, the Joule heating W^j needs to be equal or 
smaller than Wgh- 



< 



sh- 



(32) 



Therefore, we have the following constraint on the 
left hand side of Equation 



Qwhh 
Jetvp Rm 



3aP Wj 



< 



/fill /sat Wsh ' 



Jfill/sat 



/whb(/3)- (33) 



Thus, the work-heat balance poses an upper limit 
on the product of Jait/Jcqp and 1/Ru, provided 
that /fill, /sat and a are constants that do not de- 
pend on Jcrit and Rm, but only on /3. This explains 
the inverse-proportional relations observed in Fig- 
ures [3] and IH The /whb(/3) calculated with this 
interpretation using the experimental data are in 
Table 13 

We can further simplify Equation p3p by using 
the saturation predictor proposed by HGB95. The 
proposed predictors read 



aP = 0.61 ±0.06 



(34) 



S2 
877 



'16 2 

(1.21± 0.29)-27r^/— .-^0.(35) 

(36) 



Using this, Equation (|33)) is rewritten as follows: 



< 



/whb iP) 



(2.54 ±0.66) 
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1/2 



ffinfs: 



(37) 



By ignoring the dependence of fmi and /sat on /?, 
we assume fmi = 0.264 ± 0.007 and /sat = 13.1 ± 
3.1. This further simplifies the Equation p3p as: 



/whb (/?) 



(0.74 ± 0.26) /3 



1/2 



(38) 



This is in agreement with our experimental data 
(Table O within factor of 2. 

3. Distribution of The Three MRI Zones 
within The Protoplanetary Disks 

3.1. Protoplanetary Disk Model 

In the previous section, we performed the 
shearing-box simulations of MRI with nonlinear 
Ohm's law, and found the three classes of MRI 
behavior; wc named them active, dead and sus- 
tained zones. We also found the condition for the 
MRI to be sclf-sustaincd. In this section, we apply 
the findings to the global model of the protoplan- 
etary disks and analyze how they are divided into 
the three zones. 

We use Minimum- Mass Solar Nebula (M MSN) 
model introduced by iHavashi et al. I (|l985h as the 
fiducial disk model; 



S = ./sSc 



AU 



--Hi,)' 



(39) 
(40) 



Here, Sq = 1.7 x 10^ g cm^^ and To = 2.8 x 
10^ K are the surface density and the temperature 
at lAU, respectively, /s is the nondimensional 
surface density parameter. Fiducial value for the 
surface density power index is g = 3/2. Since we 
assume the isotropic equation of state (EOS), the 
ratio of specific heats 7 = 1 in our model, and the 
thermal velocities for gas molecules and electrons 
are 



Cs (r) 
Ve (r) 



(41) 
(42) 
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a 


/ \ 


/"sat 


/fill 


./whb (/?) 


400 


0.176 ±0.036 


0.222 ±0.040 


14.5 ± 1.7 


0.252 ±0.001 


7.19 ±0.81 


800 


0.143 ±0.035 


0.206 ±0.057 


17.4±2.7 


0.259 ±0.001 


9.33 ± 1.03 


1600 


0.104 ±0.027 


0.166 ±0.047 


17.4± 1.8 


0.265 ±0.002 


13.3 ±2.3 


3200 


0.0523 ±0.0206 


0.0916 ±0.0411 


12.2 ±3.9 


0.262 ±0.008 


19.0 ±3.9 


6400 


0.0236 ±0.0041 


0.0386 ±0.0072 


9.58 ± 1.04 


0.263 ±0.003 


22.2 ± 1.6 


12800 


0.0182 ±0.0032 


0.0300 ±0.0060 


9.84 ±0.88 


0.270 ±0.000 


32.6 ±2.8 


25600 


0.0185 ±0.0063 


0.0315 ±0.0128 


10.7± 1.5 


0.275 ±0.001 


58.5 ± 12.4 



Table 3: The /whb(/3) calculated from the experimental data. We first calculated the time and space 
averaged quantities a and /sat for each runs. Then for the ensemble of runs, we calculated the means and 
the standard deviations of the quantities. The ensemble constitutes of runs that (1) belong to Asustained 
zone, (2) are restarted runs {t = Wn/fl, 187r/51, 207r/f2) so that they are magnetorotationally unstable, and 
(3) have the largest product Jdit/'Ajqp ■ 1/^m so that they face the Asustained zone - xdead zone boundary. 
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Fig. 1. — Electric field amplitude as functions of 
current in the three different models of nonlinear 
Ohm's law. The symbols fid, p2, and p4 corre- 
spond to equation pT|) . and (H^]) , respectively. 



respectively, where fx is the mean molecular weight 
of the gas. 

Since we neglect the self gravity of the disk, the 
disk is in Kepler rotation and its orbital angular 
velocity is 



nir) 



(43) 



From the equilibrium between vertical pres- 
sure gradient and vertical component of the stellar 
gravity, the disk density and pressure distributions 
are: 



P (r, z) 
P{r,z) 



V2ttH 
pruH ' 



(44) 
(45) 



where H is the definition of the disk scale-height 
in this paper c.f. Equation (jl5p . 

For simplicity, we assume that the every dust 
particle to be solid sphere of the equal radius 
and density ps- The mass rud and geometrical 
cross section ad of the dust particle are 



rud 



47r , 



(46) 
(47) 
0.1 /im 



respectively. The fiducial values are ad 
and Ps = 3 g cm~'^. 

Using this rud and dust-to-gas density ratio 
fd = 0.01, the number densities of dust and gas 
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component are: 

Hg (r, z) 

nd (r, z) 



fd P 
rud 



(48) 

(49) 
(50) 



3.2. Ionization Processes 



Methods for calculating the charge equilibrium 
of the dust-pl asma in the protop lanctary disks has 
been studied (lUmebavashi k N akano 1980, 200l 
Fuiii et alj 120111) . Among those we use O09's 
method because of its numerical efficiency and 
generality. First, the gas column density above 
(Xt) and below (x^) the certain coordinate (r, z) 
in the disk are 



Xt {f, z) = 



respectively, where 
erf (x) = 



2 



2 



1 - erf 



pdz 



1 + erf 



dt 



, (51) 



, (52) 



(53) 



is the error function. 

According to SMUNOO, the effective ionization 
rate in the disk is 



CcR 

2 
exp 



exp 



Xt 
XCR 



XCR 



+ Cra, (54) 



where xcR = 96 g cm ^ CcR = 1-0 x 10 s \ 
and Cra = 6.9 x 10~23 . 

Using this, O09's nondimensional parameter 
is calculated as 



e 



(55) 



SiCsadadri^kBT' 
and r is defined as the solution of the equation 

^ 0. (56) 



1 



1 + r 



Once Equation (j56|) is numerically solved for F, 
we can calculate the number density of ions and 
electrons, ni,ne, as well as the root mean square 
of the charge per dust particle, yJJZ^e, as: 



rij (r, z) 
Tie (r, z) 



SiCsCrdUd (1 + F) ' 
Crig exp F 

SeVeCJdnd ' 

A ) ^ 2 + r A 

o2 



l + Ffld 



(57) 
(58) 
(59) 



where A 



We assume sticking probabilities Sj = 1 and = 
0.3 as assumed by O09. 

Next, we estimate the plasma conductivity us- 
ing the method of SMUNOO. The rate coeflicient 
for the collision between the neutrals and the ions 



{av), ^ 2Aln ' . (60) 

\pmH J 

We use a — 7.66 x 10~^^ cm'^ as an averaged polar- 
izability. (crw)e is the rate coefficient for collision 
between neutrals and electron, whose form is given 
in SMUNOO. The rate coefficient for dust particles 
is 



{crv)d 



An 



(61) 



This expression is valid as long as ad is much 
smaller than the mean free path of the gas 
molecules. 

With these, the magnetic diffusivity is calcu- 
lated component-wise: 



m 



c'^meng{av)e 
c^tjLmHng{av)i 
c'^fj.mHng{av)d 



(62) 
(63) 

/I 2 ' (64) 

vo = ive^'^ + vr^ + Vd^^y^ (65) 

The value of S'crit is set by the condition that the 
kinetic energy of the electrons accelerated by the 
electric field is large enough to initiate the electron 
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avalanche. 



APT/ 

E'crit 



/dp e i^'crit ^mfp; 
AW 



/dp e /nifp 
/DpeWe 



(66) 



Here, /dp — 0.43^ M'^-tr / f^e is the cocfRcicnt 
for average energy of electrons in weakly ionized 
plasma (IS05), and ^mfp — Ve/{ng{av)e) is the 
mean free path of electrons. For ionization en- 
ergy we use the value for a hydrogen molecule 
AW = 15.4 eV. 

With this, the critical current is 



J crit 



■-E'crit- 



(67) 



Note that the discharge electric field, Equation 
([66)) is calculated using the strong electric field 
limit of the electron distribution function, while 
we used the weak field limit formulae for charge 
distributions. Equations ([55|) - (|6T1) . In this paper 
we adopt this treatment for simplicity. A more 
consistent treatment will be the topic of another 
paper in preparation. 

3.3. The Self-Sustained MRI in Global 
Disk Models 

Now we study the distribution of active, sus- 
tained and dead zones in the protoplanetary disk 
models. SMUNOO gives the condition for MRI 
unstable region as follows: 



2ttv^ 

n 



< V2H /\ — < V2H. 



Combination of this with work-heat balance 
model Equation (pi]) gives the following conditions 
for active, sustained and dead zones, respectively: 



27rv^ 

n 
n 

n 



<V2Hf\^<V2H 



(69) 



27r?7 ^ /TTrr A ^crit 1 



Jcqp Rm 



</ 



whb 



(70) 



> 



V2H\/ 



2-Kl] 



> 



V2H/\ 



>/, 



whb 



(71) 



Using these Equations (|69)) - (|7T|) . we plot the 
active, sustained and dead zones for various global 
disk model Equations (pg)) - pl)) . 

First, Figure [7] shows the unstable zones for 
varying disk surface density, fs = 0.3, /s = 1 (the 
fiducial model), /e = 3, and /s = 10. In this fig- 
ure and following figures, the thick curves are the 
boundary of the work-heat balance model Equa- 
tion (|2ip , while the thin curves are the boundary of 
the instability condition in the resistive limit, i.e. 
the second condition in Equation (|68)) . The solid 
and dashed curves correspond to the plasma beta 
at the mid-plane /? = 100 and (3 ~ 1000, respec- 
tively. The active zones are marked by meshes, 
and the sustained zones are marked by horizontal 
stripes. 

Figure [S] shows the active and sustained zones 
for dust-to-gas ratio fd = fd = 0.1, fd = 0.01 
(the fiducial model), and fd ~ lO"**. The work- 
hcat balance condition Equation (PT|) is not af- 
fected by changing the dust properties such as 
dust-to-gas ratio fd or dust size a^. One can un- 
derstand this by rewriting the condition Equation 
in the following form: 



. _2 7 2 — /whb- 

47rc ^o/cqpWAz 



(72) 



This form does not include a term affected by the 
dust properties, such as the magnetic diffusivity. 
On the other hand, i^'crit is inversely proportional 
to the electron mean free path Imtp, and is propor- 
tional to the gas number density. 

In Figure [9] we study the evolution of the active 
and sustained zones as the gas density becomes 
lower while the dust density is kept constant. We 
change the set {fs,fd) from (1,0.01) (the fiducial 
model) to (0.1,0.1), (0.01,1), and (10-'^,100). In 
Figure [9] zones are marked for j3 = 1000. The mid- 
plane of the disk between the radii 2 AU — 20 AU 
becomes the sustained zone as the gas density be- 
comes 10~^ times the fiducial model. 

4. Conclusions and Discussions 

By performing numerical simulations of MHD 
with nonlinear Ohm's law of the three-dimensional 
local disks, we found hysteresis behavior for 
certain diffusivity model: If we start from the 
laminar-flow initial conditions with small seed 
fluctuations, the MRI does not activate because 
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radius (AU) radius (AU) 

Fig. 7. — Unstable regions in the protoplanetary disks. The thin sohd and thin dashed curves represent 
Xres/V^H = 1 for the cases of the magnetic field strength /? = 100 and 1000, respectively, inside of which 
is dead zone if the MRI self-sustainment is not taken into account (SMUNOO). The regions above the thick 
solid and thick dashed curves satisfies Equation ([21]) for /3 = 100 and 1000, respectively, and are sustained 
zones according to the work-heat balance model. We compare the unstable region predicted by SMUNOO to 
that predicted by our model for (3 = 1000. The unstable regions according to SMUNOO's and our model are 
marked by vertical and horizontal stripes, respectively. 
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radius (AU) radius (AU) 

Fig. 8. — Unstable regions for different dust-to-gas ratio. The MRI-unstablc region according to (SMUNOO) 
and our model are marked by vertical and horizontal stripes, respectively, for /3 — 1000. 
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radius (AU) radius (AU) 




Fig. 9. — Change of unstable regions as the gas density of the disk decreases, while the dust density is 
kept constant. The MRI-unstable region according to (SMUNOO) and our model are marked by vertical and 
horizontal stripes, respectively, for /3 = 1000. 
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of the diffusivity and the flow remains laminar; 
on the other hand, if we take MRI-turbulent 
state from an ideal-MHD simulation as initial 
conditions, MRI remains active under the same 
diffusivity model. We have surveyed in three- 
parameter space {13, Jcrit, ^^m) in search for the 
regions the self-sustained MRI in the context of 
IS05 takes place. We found the condition, the 
work-heat balance model, for this hysteresis be- 
havior to take place. The model is Wj^^ Wsh, 
where Wj is the magnetic energy dissipated by 
Joule heating per unit volume and Wsh is the 
work done by background shearing motion per 
unit volume. This leads to the proportionality 
relation j^-j^ < Uhh- 

IS05 concluded that the the energy supply from 
the shearing motion should be ~ 30 times greater 
than the energy needed to supplying the enough 
ionization for the MRI, and predicted the entire 
disk to be active. However, applying the work- 
heat balance model to various protoplanetary disk 
models, we have found that in most of the models, 
the sustained zone is above z/H > 2 — 3. 

We conclude that in the fiducial protoplane- 
tary disks environment the Joule heating (which 
has been neglected in IS05) becomes the dominant 
energy dissipation channel and constrains the self- 
sustainment of MRI, and the midplane of the disk 
remain dead. However, t he gas of the disk dis- 
sipate s ( Alexander et al. ( 2006allbl ): ISuzuki et al. 
(|2010f )) with observed timescale of 10*^ - lO'^ 



years (|Cieza et al.l (|2007t ): lHernndez et al.l (|2008t )). 
while planetesimals remain and continue planet 
formation processes. In such late phase of the 
disk, the sustained zone occupies larger volume of 
the disk. 

Although our nonlinear Ohm's law model is in- 
spired by the lightning phenomena, whether the 
Joule heating in this model takes the form of spa- 
tially and temporally concentrated stream of ion- 
izing electrons — lightning — or not, has yet to 
be studied in future works, employing nonthermal 
plasma studies. We limit ourselves to pointing 
out a few distinguishing properties of such light- 
ning which makes it an interesting subject. First, 
the work-heat balance model suggests that in sus- 
tained zones the major portion of the shearing mo- 
tion energy is converted to lightning. This means 
that the lightning is one of the most dominant en- 
ergy channel in the sustained zones. It will also 



pose a significant back-reaction to the accretion 
dynamics. We will need to reconsider the contri- 
bution of lightning in situations it has been ne- 
glected du e to lack of energy, suc h as in chondrule 
formation (j Weidenschillina 1 1 9971 ) . 

The second point is related to the redox envi- 
ronment the lightning creates. Lightning induced 
by collisional charging of water ice dust ov ercomes 
the energetics problem ( Muranushi 2010() . but if 
applied as the chondrule heating source, it suffers 
from the redox environment mismatch. Water va- 



por creates oxidizi ng environment (jClavton et al 



19811 : iRubinI 120051 ) whereas major population of 
chondrules are cons idered to have formed in re- 
ducing environment (Lofgrenll989t ConnoUv et al 



1994 bones fc DanielsonI Il997t) . However, the 
lightning in sustained zone proposed by IS05 and 
studied in this paper, is a result of pure MHD 
process, and thus is redox-neutral. Therefore, it 
can potentially explain chondrule heating in both 
reducing and oxidizing environment. 
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runID 




\ 



(1.99 ±0.60) X 10 
(1.64 ±0.69) X 10"^ 
(1.94 ±0.61) X 10^^ 
(1.30 ±0.27) X 10"^ 
(7.16 ±2.57) X 10-2 
(3.61 ±0.83) X 10-2 
(3.61 ±0.84) X 10-2 



— BxBy 



(1.05 ±0.37) X 10 
(8.27 ±3.65) X 10-2 
(9.47 ±3.19) X 10-2 
(6.08 ± 1.47) X 10-2 
(3.14 ±0.85) X 10-2 
(1.70 ±0.35) X 10-2 
(1.63 ±0.36) X 10-2 



PVxSVy 

Po 



(3.01 ± 1.49) X 10-^ 
(2.26 ±0.96) X 10-2 
(2.34 ±0.84) X 10-2 
(1.68 ±0.52) X 10-2 
(8.28 ±2.51) X 10-3 
(6.15 ±2.40) X 10 
(5.40 ± 1.84) X 10 



J2 



0.5 



(2.04 ±0.24) X 10' 
(1.94 ±0.30) X 10^ 
(2.07 ±0.26) X 10^ 
(1.85 ±0.18) X 10^ 
(1.51 ±0.15) X 10^ 
(1.24 ±0.10) X 10^ 
(1.22 ±0.09) X 10^ 



3160 
10 
220 
3180 
3612 
3616 
3620 



1-3 
,-3 



3292 





400 


1.0 


0.6 


(2.92 ± 1.56) X 10-1 


(1.67 ± 1.04) X 


10-1 


(4.81 ±3.63) X 10-2 


(2.10 ±0.43) X lOi 


1 


3293 


16T:/n 


400 


1.0 


0.6 


(2.74 ± 1.43) X 10-1 


(1.62 ±0.99) X 


10-1 


(4.63 ±3.05) X 10-2 


(2.06 ±0.48) X lOi 


(a) 


3294 




400 


1.0 


0.6 


(2.07 ±0.71) X 10-1 


(1.21 ±0.44) X 


10-1 


(3.68 ± 1.62) X 10-2 


(1.88 ±0.27) X lOi 





3295 


207r/r2 


400 


1.0 


0.6 


(2.34 ±0.82) X 10-1 


(1.36 ±0.54) X 


10-1 


(3.86 ± 1.80) X 10-2 


(2.00 ±0.30) X lOi 


J 



3352 





400 


1.0 


0.2 


(2.50 ±0.00) X 10-3 


(8.24 ±6.55) X 


10-13 


(4.83 ±2.58) X 10-^ 


(4.11 ±0.94) X 10-3 


1 


3353 




400 


1.0 


0.2 


(2.03 ± 1.27) X IQ-i 


(1.27 ±0.86) X 


10-1 


(4.07 ±2.99) X 10-2 


(1.50 ±0.60) X lOi 


(b) 


3354 




400 


1.0 


0.2 


(2.21 ±0.98) X 10-1 


(1.33 ±0.50) X 


10-1 


(4.68 ±2.68) X 10-2 


(1.66 ±0.41) X lOi 


A 


3355 


207r/r2 


400 


1.0 


0.2 


(2.52 ± 1.65) X 10-1 


(1.56 ± 1.02) X 


10-1 


(4.57 ±2.72) X 10-2 


(1.67 ±0.59) X lOi 


J 



3348 





400 


10.0 


0.2 


(2.50 ±0.00) X 10-3 


(8.24 ±6.55) X 


10-13 


(4.83 ±2.58) X 10-"^ 


(4.11 ±0.94) X 


10-3 


1 


3349 




400 


10.0 


0.2 


(2.50 ±0.00) X 10-3 


(1.77± 1.01) X 


10-^ 


(9.90 ±7.50) X 10-^ 


(3.60 ±0.54) X 


10-2 


(c) 


3350 




400 


10.0 


0.2 


(2.50 ±0.00) X 10-3 


(2.31 ±1.64) X 


10-^ 


(6.06 ±9.16) X 10-^ 


(3.49 ±0.54) X 


10-2 


X 


3351 


207r/r2 


400 


10.0 


0.2 


(2.50 ±0.00) X 10-3 


(4.85 ±3.58) X 


10-^ 


(6.62 ±8.64) X 10-^ 


(3.43 ±0.36) X 


10-2 


J 



Table 4: Statistics of the local simulations abridged. Each run is labeled by an integer. The re-start time is in the second column. Next 
three columns indicate the initial magnetic field strength, the critical current, and the magnetic Reynolds number. The physical quantity 

/ 2 2\0.5 

are represented in terms of the time average and standard deviation of the space average, i.e. A is in the format (A) ± [{A)- (A) 

. In this table are runs for ideal MHD, runs in Figure [5] that represent the behavior in (a) active, (b) sustained, and (c) dead zones, runs 
that constitute sustained-dead zone boundaries for [3 = 400, 3200 



runID 


t 


p 




i?M 






/-B.By\ 




/ pV^,SVy\ 










Jeqp 


\8ttPo/ 




\ 4^Po / 




\ Po / 






3352 





400 


1.0 


0.2 


(2.50 ±0.00) X 10" 


-3 


(8.24 ±6.55) X 10~ 


13 


(4.83 ± 2.58) X 10" 


-V 


(4.11 ±0.94) X lO"'' 


3353 




400 


1.0 


0.2 


(2.03 ± 1.27) X 10" 


-1 


(1.27 ±0.86) X 10" 


-1 


(4.07 ±2.99) X 10" 


-2 


(1.50 ±0.60) X 10^ 


3452 





400 


0.9 


0.06 


(2.50 ±0.00) X 10" 


-3 


(3.05 ±2.55) X 10"" 


13 


(4.29 ± 2.47) X 10" 


-7 


(2.22 ±0.52) X 10""^ 


3453 




400 


0.9 


0.06 


(2.50 ±0.00) X 10" 


-3 


(1.22 ±0.75) X 10" 


-8 


(2.29 ± 1.58) X 10" 


-4 


(3.02 ±0.46) X 10""2 


3448 





400 


0.3 


0.06 


(2.50 ±0.00) X 10" 


-3 


(3.05 ±2.55) X 10~ 


13 


(4.29 ± 2.47) X 10" 


-7 


(2.22 ±0.52) X 10""^ 


3449 




400 


0.3 


0.06 


(2.29 ±2.53) X 10" 


-1 


(1.43 ± 1.71) X 10" 


-1 


(4.61 ±4.95) X 10" 


-2 


(1.39 ±0.90) X 10^ 


3436 





400 


0.3 


0.02 


(2.50 ±0.00) X 10" 


-3 


(9.13 ± 12.07) X 10" 


-13 


(1.20 ± 1.54) X 10" 


-7 


(3.43 ±2.16) X 10"""^ 


3437 




400 


0.3 


0.02 


(2.50 ±0.00) X 10" 


-3 


(2.39 ± 1.08) X 10" 


~9 


(8.02 ±6.85) X 10" 


-5 


(1.17 ±0.23) X 10"2 


3432 





400 


0.1 


0.02 


(2.50 ±0.00) X 10" 


-3 


(9.55 ± 12.38) X 10" 


-13 


(8.32 ±8.82) X 10" 


-8 


(2.91 ± 1.54) X 10""^ 


3433 




400 


0.1 


0.02 


(2.00 ±2.42) X 10" 


-1 


(1.14 ± 1.50) X 10" 


-1 


(3.54 ±4.61) X 10" 


-2 


(1.23 ±0.98) X 10^ 


3372 





3200 


1.0 


0.2 


(1.11 ±0.62) X 10" 


-1 


(5.22 ±2.79) X 10" 


-2 


(1.39 ±0.68) X 10" 


-2 


(1.58 ±0.31) X 10^ 


3373 


167r/r2 


3200 


1.0 


0.2 


(1.00 ±0.41) X 10" 


-1 


(4.29 ± 1.61) X 10" 


-2 


(9.90 ±3.04) X 10" 


-3 


(1.52 ±0.24) X 10^ 


3460 





3200 


0.9 


0.06 


(3.12 ±0.00) X 10" 


-4 


(3.02 ± 1.86) X 10" 


12 


(4.68 ±2.21) X 10" 


-7 


(1.89 ±0.39) X 10"^ 


3461 


167r/r2 


3200 


0.9 


0.06 


(3.35 ± 0.24) X 10" 


-4 


(9.82 ± 10.64) X 10 


-6 


(3.36 ±3.24) X 10" 


-5 


(5.07 ±2.78) X 10"2 


3456 





3200 


0.3 


0.06 


(3.12 ±0.00) X 10" 


-4 


(3.02 ± 1.86) X 10" 


12 


(4.68 ±2.21) X 10" 


-7 


(1.89 ±0.39) X 10"^ 


3457 




3200 


0.3 


0.06 


(1.19 ±0.48) X 10" 


-1 


(5.44 ± 1.88) X 10" 


-2 


(1.42 ± 0.51) X 10" 


-2 


(1.65 ±0.26) X 10^ 


3444 





3200 


0.3 


0.02 


(3.12 ±0.00) X 10" 


-4 


(8.60 ±6.12) X 10" 


13 


(5.18 ±2.54) X 10" 


-7 


(1.55 ±0.35) X 10"^ 


3445 




3200 


0.3 


0.02 


(9.38 ±4.06) X 10" 


-2 


(4.50 ± 1.88) X 10" 


-2 


(1.15 ±0.46) X 10" 


-2 


(1.29 ±0.27) X 10^ 


3440 





3200 


0.1 


0.02 


(3.12 ±0.00) X 10" 


-4 


(8.60 ±6.12) X 10" 


13 


(5.18 ±2.54) X 10" 


-7 


(1.55 ±0.35) X 10"^ 


3441 




3200 


0.1 


0.02 


(8.15 ±3.91) X 10" 


-2 


(3.79 ± 1.70) X 10" 


-2 


(1.05 ±0.50) X 10" 


-2 


(1.41 ±0.23) X 10^ 



Table 4: (continued) 



